      subroutine propgeom(x,y,N,NFIX,compr,area,xc,yc,Ju,Jv,Juv)

      include 'param_dim.inc'

      integer N,NFIX
      double precision x(N),y(N),area,area_parc,gxi,gyi
      double precision gx_sum,gy_sum,compr
      double precision Jx,Jy,Jx_parc,Jy_parc,Jxy,Jxy_parc,Ju,Jv,Juv
      double precision xc,yc,xc1,yc1,xi,yi,xf,yf,alfa

    
      area=0
      compr=0
      gx_sum=0
      gy_sum=0
      Jx = 0
      Jy = 0
      Jxy = 0
      
      do i=NFIX+1,N-1
         xi=x(i)
         yi=y(i)
         xf=x(i+1)
         yf=y(i+1)

         area = area + area_parc(xi,yi,xf,yf)
         compr= compr + sqrt((xf-xi)*(xf-xi)+(yf-yi)*(yf-yi))
         gx_sum = gx_sum + gxi(xi,yi,xf,yf)
         gy_sum = gy_sum + gyi(xi,yi,xf,yf)

      enddo

      xi=x(N)
      yi=y(N)
      xf=x(NFIX+1)
      yf=y(NFIX+1)

      area =-(area + area_parc(xi,yi,xf,yf))
      compr= compr + sqrt((xf-xi)*(xf-xi)+(yf-yi)*(yf-yi))

      gx_sum = -(gx_sum + gxi(xi,yi,xf,yf))
      gy_sum = -(gy_sum + gyi(xi,yi,xf,yf))

      xc = gx_sum/area
      yc = gy_sum/area
      

      
      do i=NFIX+1,N-1
         xi=x(i)  -xc1
         yi=y(i)  -yc1
         xf=x(i+1)-xc1
         yf=y(i+1)-yc1
         
         Jx = Jx + Jx_parc(xi,yi,xf,yf)
         Jy = Jy + Jy_parc(xi,yi,xf,yf)
         Jxy= Jxy+ Jxy_parc(xi,yi,xf,yf)
      enddo 
      
      xi=x(N)    -xc1
      yi=y(N)    -yc1
      xf=x(NFIX+1)-xc1
      yf=y(NFIX+1)-yc1
      
      Jx = -(Jx + Jx_parc(xi,yi,xf,yf))
      Jy = -(Jy + Jy_parc(xi,yi,xf,yf))  
      Jxy=-(Jxy+ Jxy_parc(xi,yi,xf,yf))

c     Calculo dos momentos de inercia segundo as direcoes...
c     principais 
c     lembrar que alfa=2*teta
c     Juv so para confirmar que da zero


        alfa=atan(-2*Jxy/(Jx-Jy))
	if(Jx .eq. Jy )then 
		alfa=0.
	endif


      Ju = (Jx+Jy)/2 + (Jx-Jy)*cos(alfa)/2 - Jxy*sin(alfa)
      Jv = (Jx+Jy)/2 - (Jx-Jy)*cos(alfa)/2 + Jxy*sin(alfa)
      Juv= (Jx-Jy)*sin(alfa)/2 + Jxy*cos(alfa)


      return
      end

      function area_parc(xi,yi,xf,yf)
      double precision area_parc,xi,yi,xf,yf
      
      area_parc=(yf* xi - xf*yi)/2
      
      return
      end

      
	function gxi(xi,yi,xf,yf)
	double precision gxi,xi,yi,xf,yf
	
		gxi=(xf*yf*xi - xf*xi*yi)/6 - (xf*xf*yf)/3
     &             +(yf*xi*xi)/6 - (xf*xf*yi)/6 + (xi*xi*yi)/3
	
	return
	end

	function gyi(xi,yi,xf,yf)
	double precision gyi,xi,yi,xf,yf
	
		gyi=(yf*xi*yi - xf*yf*yi)/6 + (xf*yf*yf)/3 - (xf*yi*yi)/6 
     &            + (yf*yf*xi)/6 - (xi*yi*yi)/3
	
	return
	end

	function Jx_parc(xi,yi,xf,yf)
	double precision Jx_parc,xi,yi,xf,yf,xc,yc
	
	Jx_parc= (xf*yf*yf*yf)/4 - (xf*yi*yi*yi)/12 + (yf*yf*yf*xi)/12
     &         - (xi*yi*yi*yi)/4 - (xf*yf*yi*yi)/12 - (xf*yf*yf*yi)/12
     &         + (yf*xi*yi*yi)/12 + (yf*yf*xi*yi)/12
	
	
	return
	end

      function Jy_parc(xi,yi,xf,yf)
      double precision Jy_parc,xi,yi,xf,yf,xc,yc
    
      Jy_parc= (yf*xi*xi*xi)/12 - (xf*xf*xf*yf)/4 - (xf*xf*xf*yi)/12 
     &       + (xi*xi*xi*yi)/4 + (xf*yf*xi*xi)/12 + (xf*xf*yf*xi)/12
     &       - (xf*xi*xi*yi)/12 - (xf*xf*xi*yi)/12
      
      return
      end

      function Jxy_parc(xi,yi,xf,yf)
      double precision Jxy_parc,xi,yi,xf,yf

      Jxy_parc= (xf*yf*yf*xi)/12 - (xf*xf*yf*yi)/12 - (xf*xi*yi*yi)/12 
     &        + (yf*xi*xi*yi)/12 + (xf*xf*yf*yf)/8  - (xf*xf*yi*yi)/24 
     &        + (yf*yf*xi*xi)/24 - (xi*xi*yi*yi)/8
      
      return
      end
      
      subroutine centroide(x,y,N,NFIX,xc,yc,area)

      include 'param_dim.inc'

      integer N,NFIX
      double precision x(NX),y(NX),area,area_parc,gxi,gyi,gx_sum
      double precision xc,yc,xi,yi,xf,yf
      
      area=0
      gx_sum=0
      gy_sum=0

      
      do i=NFIX+1,N-1
         xi=x(i)
         yi=y(i)
         xf=x(i+1)
         yf=y(i+1)

         area = area + area_parc(xi,yi,xf,yf)
         gx_sum = gx_sum + gxi(xi,yi,xf,yf)
         gy_sum = gy_sum + gyi(xi,yi,xf,yf)
         
      enddo

      xi=x(N)
      yi=y(N)
      xf=x(NFIX+1)
      yf=y(NFIX+1)

      area =(area + area_parc(xi,yi,xf,yf))
      gx_sum = gx_sum + gxi(xi,yi,xf,yf)
      gy_sum = gy_sum + gyi(xi,yi,xf,yf)
     
      xc = gx_sum/area
      yc = gy_sum/area
     
      
      return
      end

